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ABSTRACT 

The theory of shock acceleration predicts the maximum particle energy to be 
limited only by the acceleration time and the size (geometry) of the shock. This 
led to optimistic estimates for the galactic cosmic ray energy achievable in the 
SNR shocks. The estimates imply that the accelerated particles, while making 
no strong impact on the shock structure (test particle approach) are neverthe¬ 
less scattered by strong self-generated Alfven waves (turbulent boost) needed to 
accelerate them quickly. We demonstrate that these two assumptions are in con¬ 
flict when applied to SNRs of the age required for cosmic ray acceleration to the 
“knee” energy. 

We study the combined effect of acceleration nonlinearity (shock modification 
by accelerated particles) and wave generation on the acceleration process. We 
show that the refraction of self-generated waves resulting from the deceleration 
of the plasma flow by the pressure of energetic particles causes enhanced losses 
of these particles. This effect slows down the acceleration and changes the shape 
of particle spectrum near the cut-off. The implications for observations of TeV 
emission from SNR remnants are also discussed. 

Subject headings: acceleration of particles—cosmic rays —shock waves—supernova 
remnants—turbulence 


1. Introduction 

The first-order Fermi or diffusive shock acceleration (DSA) has been long considered as 
responsible for the production of galactic cosmic rays (CRs) in supernova remnants (SNRs), 
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as well as for the radio, x- and 7 -ray emission from these and other shock related objects. 
The most crucial characteristic of this process that is usually examined in terms of its 
capability to explain a given observation, is the rate at which it operates. Indeed, what 
is often expected from the theory or even inferred from the observations is an extended 
particle energy spectrum, frequently a power-law, but more rapidly decaying at the highest 
energies observed. Often, this decay is referred to as an energy or momentum cut-off and 
is usually associated with the finite acceleration time or with losses if their rate exceeds the 
acceleration rate. As long as the losses are unimportant, the cut-off p ma x(t) advances with 
time according to the following equation 
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whereas in the presence of losses the acceleration rate Pmax/tacc may be equated to the loss 
rate to yield a steady state value of p max ■ The acceleration time scale is determined by (e.g., 
Axford 1981) 
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with U\ and u 2 being the upstream and downstream flow speeds in the shock frame and with 
K\ )2 being the particle diffusivities in the respective media. One may recognize in the last 
formula the sum of average residence times of a particle spent upstream and downstream 
of the shock front before it completes one acceleration cycle, integrated over the entire 
acceleration history from p m in to p m ax■ Given the flow speeds u\ )2 which, in many cases 
are known reasonably well, the most sensitive quantity is the particle diffusivity k. This, in 
turn, is determined by the rate at which particles are pitch angle scattered by the Alfven 
turbulence. If the latter was just a background turbulence in the interstellar medium (ISM), 
the acceleration process would be too slow to produce the galactic CRs in SNRs (e.g., Lagage 
& Cesarsky 1983). However it was realized (e.g., Bell 1978; Blandford & Ostriker 1978) 
that accelerated particles should create the scattering environment by themselves generating 
Alfven waves on the cyclotron resonance kpn/m = cn C j, where k is the magnitude of the 
wave vector (directed along the magnetic field), p, //, m and c o ci are the particle momentum, 
the cosine of its pitch angle, mass and non-relativistic (eB/me) gyro-frequency. Note, that 
the diffusive character of particle transport (and determination of k) has been rigorously 
obtained within a quasi-linear theory, i.e., it is subject to constraints on the turbulence 
level. 


The wave generation, however, proved to be very efficient (see e.g., Volk et al. 1984 and 
the next section). In particular, using, again, the quasi-linear approximation, the normal- 
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ized wave energy density ( SB/Bo) 2 may be related to the partial pressure P c of CRs that 
resonantly drive these waves through 

(SB/B 0 f ~ M A P c /pu 2 (3) 

where M A is the Alfven Mach number and pu 2 is the shock ram pressure. Since M A is 
typically a large parameter, SB/Bo may become larger than unity even if the acceleration 
itself is relatively inefficient, i.e., if P c /pu 2 <C 1. Strictly speaking this invalidates the quasi- 
linear approach as a means for describing the generation of strong turbulence at shocks. 
The commonly accepted way to circumvent this difficulty is to assume that the turbulence 
saturates at SB/Bo ~ 1, which means that the m.f.p. of pitch angle scattered particles is 
of the order of their gyro-radius r g . Then, k = n B = cr g (p)/ 3, where the speed of light 
c is substituted instead of CR velocity and kb stands for the Bohm diffusion coefficient. 
This immediately sets the acceleration time scale (2) at the level of particle gyro-period 
(i eB/p ) 1 times ( c/u \) 2 . In principle, the turbulence level SB/Bo significantly exceeding 
unity is possible in local shock environments (see e.g., numerical studies by Bennett & Ellison 
1995 and Bell & Lucek 2000). As a consequence of that the diffusion coefficient could be 
even smaller than kb, and hence, the acceleration rate would be faster than it is commonly 
believed to be. At the same time, since usually Alfvenic type turbulence is considered, 
the respective velocity perturbations must be super-Alfvenic and supersonic, which raises 
questions about its ability to sustain itself in an extended area without rapid dissipation 
that will decrease the SB/B 0 level. Likewise, decreasing of turbulence level below the Bohm 
limit, for example due to the finite extent of the turbulence zone upstream, should slow down 
the acceleration (Lagage & Cesarsky 1983). 

However the acceleration rate given by eq.(2) with k = kb was found to be fast enough 
to explain (at least marginally) the acceleration of CRs in SNRs up to the “knee” energy 
~ 10 15 eE over their life time. Much further optimism has been caused by the studies 
of Drury et al. (1994) and Naito & Takahara (1994). They analyzed the prospects for 
detection of super-TeV emission from nearby SNR that should be produced by the decays 
of 7T° mesons born in collisions of shock accelerated protons with the nuclei of interstellar 
gas. The expected fluxes were shown to be detectable by the imaging Cerenkov telescopes. 
Moreover, the EGRET (Esposito et al. 1996) detected a lower energy (< GeV) emission 
coinciding with some galactic SNRs. The spectra also seemed consistent with the DSA 
predictions. One may even argue that the low energy EGRET data verified one of the most 
difficult elements of the entire acceleration mechanism, the so called injection. In essence, 
this is a selection process (not completely understood) whereby a small number of thermal 
particles become subject to further acceleration (see Gieseler et al. 2000; Zank et al. 2001 for 
the latest development of the injection theory and Maikov & Drury 2001 for a review) and 
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may be then treated by standard means of the DSA theory that was designed to describe 
particles with velocities much higher than the shock velocity. Therefore, what seemed left 
for the theory was to continue the EGRET spectrum (that sets the normalization constant, 
or injection rate) with some standard DSA slope (nearly E~ 2 or somewhat steeper) and to 
predict the 7 -ray flux in the TeV range where it could be detected by the Cerenkov telescopes. 

Unfortunately, despite the physical robustness of the arguments given by Drury et al. 
(1994); Naito & Takahara (1994), no statistically significant signal that could be attributed 
to any of the EGRET sources was detected. The further complication is that some critical 
energy band between GeV and TeV energies is currently uncovered by available instruments. 
Therefore, based on these observational results it was suggested (e.g., Buckley et al. 1998) 
that there is probably a spectral break or even cutoff somewhere within this band. However 
the spectrum above GeV energies remains an enigma. This will be resolved perhaps with 
the launch of the GLAST mission and when the new generation of Cerenkov telescopes with 
lower energy thresholds begin to operate. However, the discovery of the 100 TeV emission 
from SNR 1006 (Tanimori 1998), as well as some other remnants not seen by the EGRET 
at lower energies (see, e.g., Aharonian et al. 2001; Allen et al. 2001; Kirk & Dendy 2001 
for a complete discussion), although almost universally identified with electrons diffusively 
accelerated to similar energies, is widely interpreted as a strong support of the mechanism 
itself. The above suggests, however, that in reality it might be not as robust as is its simplified 
test particle version with enhanced turbulence and particle scattering. 

In this paper we attempt to understand what may happen to the spectrum provided 
that the acceleration is indeed fast enough to access the TeV energies over the life time 
of SNRs in question. Our starting point is that the fast acceleration also means that the 
pressure of accelerated particles becomes significant in an early stage of Supernova evolution 
so that the shock structure is highly nonlinear. At the first glance this should not slow down 
acceleration since, according to eq.(3), this changes the turbulence level thus improving 
particle confinement near the shock front and thus making acceleration faster (smaller k). 
However, the formation of a long CR precursor (in which the upstream flow is gradually 
decelerated by the pressure of CRs, P c ) influences the spectral propeHies of the turbulence 
by affecting the propagation and excitation of the Alfven waves. This effect is twofold. First 
the waves are compressed in the converging plasma flow upstream and are thus blue-shifted, 
eliminating the long waves needed to keep exactly the highest energy particles diffusively 
bound to the accelerator. Second, and as a result of the first, at highest energies there 
remain fewer particles than expected so that the level of resonant waves is smaller and hence 
the acceleration rate is lower. We believe that these effects have been largely overlooked 
before which may have substantially overestimated the particle maximum energy in strongly 
nonlinear regimes. 
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2. Basic Equations and Approximations 

We use the standard diffusion-convection equation for describing the transport of high 
energy particles (CRs) near a CR modified shock. First, we normalize the distribution 
function f(p) to p 2 dp\ 


, v dl _ d_ K dl = 19U W 

dt dx dx K dx 3 dx P dp 

Here x is directed along the shock normal which for simplicity is assumed to be the direction 
of the ambient magnetic held. The two quantities that control the acceleration process are 
the how prohle U(x) and the particle diffusivity hi(x,p). The hrst one is coupled to the 
particle distribution / through the equations of mass and momentum conservation 
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is the pressure of the CR gas, P g is the thermal gas pressure, p is its density. The lower 
boundary in the momentum space, pi n] separates CRs from the thermal plasma that is not 
directly involved in this formalism but rather enters the equations through the magnitude 
of / at p = pi n j which specihes the injection rate of thermal plasma into the acceleration 
process. The particle momentum p is normalized to me. 


Since we shall be primarily concerned with the wave generation and particle confinement 
upstream of the discontinuity, we assume that the upstream region is at x > 0, so that the 
velocity prohle can be represented in the shock frame as U(x) = —u(x) where the (positive) 
how speed u(x) jumps from u 2 = w(0—) downstream to u 0 = u(0+) > u 2 across the sub-shock 
and then gradually increases up to Ui = u(+ oo) > u 0 (see Figure 2a below). 

We may neglect the contribution of the gas pressure to eq.(6) in the upstream region 
(x > 0, but not at x < 0) restricting our consideration to the high Mach number shocks, 
M 3> 1. The gas pressure is retained when treating the sub-shock (discontinuous part of 
the shock structure) which, however, can be simply described by the conventional Rankine- 
Hugoniot jump condition (since P c does not vary on this scale) 



6 


7+1 , ( 8 ) 

it 2 7 - 1 + 2 Mq 2 

Here Mq is the Mach number in front of the sub-shock. When the flow compression in the 
CR precursor can be considered as adiabatic, this can be expressed through the given far 
upstream Mach number in a standard way, Mg = M 2 /R 1+1 , where R = ui/uo is the flow 
precompression in the CR precursor. We shall also set 7 = 5/3 in what follows. 

Turning to determination of the CR diffusion coefficient k, we note that since the CR 
precursor scale height is ~ n(p m ax)/u 1 ~ (c/ui)r g (p max ), which is cju\ 1 much larger 
than the longest wave in the spectrum ~ r g (p max ) we can use a wave kinetic equation in the 
eikonal approximation for describing the evolution of Alfven waves 


dN k du dN k 
dt dk dx 


du d N k 
dx dk 


lk N k + St{N k } 


( 9 ) 


Here N k is the number of wave quanta and ix is the wave frequency x = —ku + kV a — —ku. 
The left hand side has a usual Hamiltonian form that states the conservation of N k along the 
lines of constant frequency ix(k,x) = const on the k,x plane. The first term on the r.h.s. 
describes the wave generation on the cyclotron instability of a slightly anisotropic particle 
distribution. It can be expressed through its spatial gradient. The resonance condition 
for the wave-particle interaction contains also the particle pitch angle cos ^ 1 p by means of 
the following expression kpp = eB/c which, generally speaking, requires the treatment of 
particle distribution in two dimensional momentum space (p, p). A significant simplification 
can be achieved by the so called “resonance sharpening” procedure (Skilling 1975; Drury et 
al. 1996) whereby a certain “optimal” value of p is ascribed to all particles and the resonance 
condition puts k and p into a one-to-one relation, i.e., kp = const . The second term on 
the r.h.s. stands for nonlinear wave-particle and wave-wave interactions such as the induced 
scattering of waves on thermal protons and mode coupling (Sagdeev & Galeev 1969). We 
will suggest a simple model for this nonlinear term in section 3.2. 


To conclude this subsection we emphasize that while eqs.(4-6) already treat the acceler¬ 
ation process and flow structure on equal footing, the fluctuation part given by eq.(9) must 
be included in this treatment and, as we shall see in the sequel, it by no means plays a 
subdominant role in this triad. 
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2.1. The significance of acceleration nonlinearities 

There are two aspects of the acceleration where nonlinearity is crucial for its outcome. 
The first aspect is the excitation of scattering waves by accelerated particles and the second 
one is the backreaction of these particles on the shock structure. The latter is critical both 
for the particle injection and wave excitation that is, to particle confinement. 

Indeed, the system 4-7 self-consistently describes particle acceleration and the shock 
structure (nonlinearly modified by the particle pressure) only if the particle scattering law 
is known (which is contained in the diffusion coefficient k) and the injection rate from the 
thermal plasma is also known (the normalization of the particle distribution f(p) in eq.(4)). 
Physically, the scattering rate determines the particle maximum momentum p ma x , as eq.(2) 
indicates. The difficulty, however, is that both the cut-off momentum p m ax and the wavenum¬ 
ber cutoff of the scattering turbulence change in time simultaneously (one controlling the 
other) due to the cyclotron resonance condition. However the speed at which they change 
has not been calculated self-consistently. The linear solution given by eqs.(l) and (2) is 
essentially based on the assumption that p m ax is growing due to already existing stationary 
turbulence. In reality, the particle energy cut-off and the corresponding cut-off on the wave 
spectrum, as we mentioned, both advance together and, since waves need to grow from a 
very small background amplitude at each current cut-off position, an additional slow down 
must be introduced in the entire process. A good analogy here is the problem of beam 
relaxation in plasmas (Ivanov 1978) (where a front on the particle velocity distribution also 
propagates on self-generated rather than on pre-existing resonant waves). This suggests that 
the speed of the front in momentum space, as given by eqs.(l) and (2) should be reduced 
by a factor ~ \yi{W/Wism), where H'Vsm is the background turbulence amplitude and W 
is the saturated wave amplitude generated by accelerated particles. As we mentioned, the 
latter may be associated with Wb = Bq/8tt so that the acceleration time given by eq.(2) 
may increase by a factor ~ 10 (e.g., Achterberg et al. 1994 estimate Wism/XVb ~ ICE 5 ). 
Evidently, the additional In -factor takes care of the time needed for waves to grow before 
they start to scatter particles with current momentum p at the Bohm rate. 

The above consideration also shows that particles with p < p max are confined to the 
shock through fast pitch angle scattering while particles with p > p max are scattered only very 
slowly due to the absence of self-generated waves and leave the accelerator. Mathematically, 
this means f(p > p ma x) = 0 or n(p > p max ) = oo. Note that the propagating front solution 
must produce different (sharper) cut-off shape at p = p max (t) than approaches based on the 
pre-existing turbulence, i.e., on a prescribed (for all p) n(p), (e.g., Berezhko et al. 1996). 
Even if the speed and the form of the front at p ma x{t) are unknown, the above ansatz allows 
analytic solution of the system (4-7) (Maikov 1997) for p < p max in the limit of strong 



shocks M 3> 1, high maximum momentum p max (that may slowly advance in time) and for 
essentially arbitrary, in particular, Bohm n(p) dependence for p < p max (as we mentioned, 
often assumed in numerical studies, however, for all p, e.g., Duffy 1992). The analytic 
solutions are tabulated e.g., by Maikov & Drury (2001) and extensively used below. 

Since waves are generated by accelerated particles upstream in the precursor, the main 
nonlinear impact on the wave dynamics and thus on p max must be from the flow pre¬ 
compression. The latter can be characterized by the parameter R = U\/uq which is shown in 
Figure 1 as a function of injection parameter v for different maximum momenta p max - The 
injection parameter v is related to the normalization of particle distribution function / in 
eq.(4) as 


v = 


47 t me 2 


ptnjMP 


mj) 


( 10 ) 


where fo(p) is the downstream value of /. In this form the injection rate v naturally appears 
as a coefficient in front of the CR pressure in the momentum flux conservation eq.(6) when 
the CR pressure is normalized to the ram pressure p\u\ (see eq.[12] in the next section). 


One aspect of the solution shown in Figure 1 which is important here is that for any 
given injection rate u, the growing maximum momentum p max (t ) will ultimately exceed a 
critical value, beyond which the test particle regime fails to exist. (It is natural to assume 
that the acceleration starts at this regime i.e., where R « 1, e.g., point A on Figure 1). 
Formally, the system must then transit to a much higher R that will be still very sensitive 
to the current values of i/(fa 0.01) and p max (— 10 6 ) as may be seen from Figure 1 (point 

B) . Obviously, the further development of the acceleration process will depend on how the 
parameters v and p max react to this strong increase of R. One possibility is to assume that 
simply a constant fraction of the sub-shock plasma is injected so that the injection rate 
substantially increases because the plasma density at the sub-shock grows linearly with R. 
Then, the system must leave the critical region where the R{v) dependence is very sharp or 
even nonunique and proceed to a highly supercritical regime characterized by higher v (point 

C) . The curve R[y) saturates there at the level oc M 3//4 which in the most straightforward 
way may be deduced from the condition of the sub-shock preservation, M 0 > 1, (see eq.[8]). 
A general formula for R(u, M) with the M 3 / 4 scaling as a limiting case may be found in 
(Maikov 1997). This scenario was realized in many numerical models (e.g., Ellison & Eichler 
1985; Kazanas & Ellison 1986; Berezhko et al. 1996), since they normalized the injection 
parameter to the plasma density at the sub-shock p 0 = piR, which should clearly lead to the 
R oc M 3 / 4 scaling. Obviously, the pre-compression R and thus the acceleration efficiency will 
then be insensitive to v (in deep contrast to the case v ~= const , point B) since the point 
C is already on the saturated part of the R{v) curve. Often, this insensitivity is observed 
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in numerical studies with the parameterized injection rate (e.g., Berezhko et ah 1996), so it 
is tempted to conclude that we do not need to know the injection rate very accurately, as 
soon as it exceeds the critical rate. 

However, the injection rate is known to be suppressed by a number of self-regulating 
mechanisms such as trapping of thermal particles downstream by the injection driven tur¬ 
bulence (Maikov 1998) and cooling of the downstream plasma in strongly modified shocks. 
These effects are believed to more than compensate the compressive growth of plasma density. 
Recently, these effects have been systematically included in numerical studies by Gieseler et 
al. (2000); Kang et al. (2001). They did not confirm the simple v oc R rule. Instead they 
indicate that in course of nonlinear shock modification accompanied by growing R , the in¬ 
jection rate v remains remarkably constant (Gieseler et al. 2000). Moreover, the preliminary 
results of a new adaptive mesh refinement (AMR) modification of these scheme allowing 
higher p m axj indicate that the injection efficiency may even begin to decrease with growing 
Pmax (Kang et al. 2001, 2002). These self-regulation mechanisms are applicable to both 
strictly parallel and oblique shocks of which the former ones is clearly an exceptional case. 
Even slightly oblique shocks have an additional self-regulation of injection via a nonlinearly 
increasing obliquity. Indeed, since the tangential magnetic field component B t is amplified 
at the sub-shock by the factor of R , the sub-shock may be strongly oblique even if the shock 
itself is not. This leads to exponentially strong suppression of leakage of downstream thermal 
particles upstream (for a Maxwellian downstream distribution) since the intersection point 
of a held line (the particles sit on) with the shock front rapidly escapes from these particles. 
On the other hand, enhanced particle reflection off the oblique subshock should increase 
injection. 

The inspection of Figure 1 shows that if we (conservatively) assume n(R) = const (AB) 
rather than u oc R (AC) the results will differ dramatically particularly in terms of the 
injection rate. Note that particle spectra that correspond to the points A and C also differs 
very strongly (see Maikov & Drury 2001 for graphical examples). What is important for the 
subject of the present paper is that in both these cases as well as for any other point on the 
part BC of the bifurcation curve, the compression R is very high. It has been pointed out 
by Maikov et al. (2000), that this must have strong impact not only on the injection rate 
as discussed above, but also on the wave propagation and thus on the particle confinement. 
This in turn should lead to significant reduction of the maximum momentum achievable by 
this acceleration mechanism. We will quantify these effects in the next section. 
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3. Analysis 

Returning to eqs.(4) and (9). it is convenient to use the wave energy density normalized 
to d In k and to the energy density of the background magnetic field Bq/8tt instead of Nk 


h = 


k 2 V A 
Bit 8* 


N, 


( 11 ) 


along with the partial pressure of CRs normalized to d In p and to the shock ram pressure 


piui 


p 


47 r me 2 p 5 

Ttw I ^rt f(p,x 
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Using these variables, denoting g = P/p, assuming a steady state and p> 1, eqs.(4,9) can 
be rewritten as (Bell 1978; Drury et ah 1996) 


d ( dg \ 1 dg 

— [ u g + k— = -u x p— 

ox \ ox J 3 op 

(13) 

dI o d I 2 u\ d 

U dx +UxP dpp 2 ~ V A dx P ^ ^ 

(14) 


Here u x = du/dx and the wave intensity I = I(p) = Ik is now treated as a funcion of p 
rather than k according to the resonance relation kp = const. The CR diffusion coefficient 
k can be expressed through the wave intensity by 



(15) 


where Kr ip) is the Bohm diffusion coefficient. The difference between these equations and 
those used by, e.g., Bell (1978); Drury et ah (1996) is due to the terms with u x ^ 0 and the 
»57-term on the r.h.s. of eq.(14). Far away from the sub-shock where u x —> 0, and where the 
wave collision term is also small due to the low particle pressure P, one simply obtains 


I 


2ui 

V A 


(16) 


Note, that this shows the limitation of the linear approach in the case of strong shocks 
Ma = Ui/Va 3> 1. The most important change to the acceleration process comes from the 
terms with u x ^ 0. Indeed, let us recall first how the equation (13) may be treated in the 
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linear case u x = 0 for x > 0. We integrate both sides between some x > 0 and +oo, which 
yields 


uig + 


kqVa i dg 

2 u\ g dx 


0 


( 17 ) 


where we denoted = kb/p — const for 1. Although this equation has a formal spatial 
scale l ~ Ko/uiMAg, its only solution is a power law 


9 oc 1 / (x + x 0 ) (18) 

and thus has no scale (xo = xo(p) is an integration constant). It simply states the balance 
between the diffusive flux of particles upstream (second term in eq.[17]) and their advection 
with thermal plasma downstream (the first term). As we shall see, this balance is possible 
not everywhere upstream and the physical reason why it appears to be so robust in the case 
u x = 0 is that flows of particles and waves on the x, p-plane (including the diffusive particle 
transport) are both directed along the x-axis. If, however, the flow modification upstream 
is significant (u x > 0, x > 0), the situation changes fundamentally. Figure 2 explains how 
the flows of particles and waves on the x, p-plane become misaligned even though they are 
both advected with the thermal plasma. In fact, the flows separate from each other and, 
since neither of them can exist without the other (waves are generated by particles that, in 
turn, are trapped in the shock precursor by the waves) they both disappear in some part of 
the phase space. To understand how this happens we rewrite eqs. (13-14) in the following 
characteristic form (we return to the particle number density /) 

( d 1 d \ . d d f 

d \ I 2 u] d 1 , , 

+ UxP dp)f = v^~d~x P ~ f St { } 

One sees from the l.h.s.’s of these equations that particles are transported towards the sub¬ 
shock in x and upwards in p along the family of characteristics up 3 = const, whereas waves 
move also towards the sub-shock but downwards in p along the characteristics u/p = const. 
As long as u(x) does not significantly change the waves and particles propagate together 
(along x-axis) as e.g., in the case of unmodified shock or far away from the sub-shock where 
u x —> 0. When the flow compression becomes important (u x ^ 0) their separation leads to 
decrease of both the particle and wave energy densities towards the sub-shock. To describe 
this mathematically, let us assume first that the relation (16) between P and / is still a 



(19) 

( 20 ) 
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reasonable approximation even if u x is nonzero but small. Then, integrating eq.(13) again 
between some x > 0 and x — oo, instead of (17) we obtain 


ug + u i 


(7 (9a: 


1 

3 



dg 

u x p—dx 

op 


( 21 ) 


In contrast to the solution of eq.(17) the length scale L = k 0 /2uiMa enters the solution of 
this equation. This is because it has a nonzero r.h.s. In the next subsection we obtain the 
solution of this equation that rapidly changes on a scale ~ L. 


3.1. Internal asymptotic solution for g 

First we note that L < f c where L c = K,(p max )/ui is the total scale height of the CR 
precursor on which u(x) changes. Next, in addition to x and p, we introduce a fast (internal) 
variable £(x,p) as follows 


£ 


X — Xf (p) 

L 


( 22 ) 


where x = Xf(p) is some special curve on the x,p plane which bounds the solution and will 
be specified later. We rewrite eq.(21) for 


£ — fixed, L —> 0 

Separating fast variable terms on the r.h.s. by replacing 


(23) 


dg/dp —> dg/dp — L 1 (dxf/dp) (dg/d £), 
to the leading order in L/L c —> 0, we obtain 


Mi dg 1 dip 1 f°° dG 

Mf q H-— = -v—— (Ct — q) -/ u x p——dx 

g 3 y dp K 3 J x dp 


(24) 


Here we denoted Uf(p) = u [xf(p)] and 


G{x,p) = lirn g{i,x,p) (25) 

KX) 

The existence of this limit will be confirmed upon obtaining the solution of eq.(24) below. 
First, we introduce the following notations 
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, . 1 du{ 

w[p) =u { + -p— 
3 dp 

, 1 du t 1 /•“ 

S{X ’ P) = 3 P 1^ G ~ 3 L 


dG , 

u x p—dx 

op 


Eq.(24) then can be rewritten as: 


wg 


Mi dp _ Q 
9 ' 


and its solution can be thus written as 


(26) 


9{£,x,p) 


S(x,p ) 

w(p) + e~ s ^/ ui 


(27) 


One sees that the limit in eq.(25) indeed exists and is equal to G — S/w. Furthermore, 
eq.(27) describes a transition front on the particle distribution between its asymptotic value 
g = G at £ —> oo and g = 0 at £ —> —oo. This front solution establishes as a result of particle 
losses caused by the lack of resonant waves towards the sub-shock as we argued discussing 
eqs.(19,20). Note that according to the ordering in eq.(23) we should set x = X{(p) in S(x,p ) 
when solving (26) for gr(£) and we indeed must do it for £ ~ 1 as well as for all negative 
£ < 0. In the limit £ —> oo, however, we may use the result (27) for arbitrary x > Xf(p) 
since it remains valid in this case, however, it merely states that in this region the complete 
solution is represented by its “external” part G(x,p ) (eq. [25]). This, in turn, is yet to be 
determined. Before we do this in section 3.3 we should verify the validity of the internal 
solution and calculate its unknown function Xf(p). 


3.2. Nonlinear modification of the internal solution. Determination of Xf (p) 

The way we resolved eq.(14) for I (see eq.(16)) may become inadequate for two reasons. 
First, the second term on the l.h.s. of eq.(14) may become comparable with the first one. 
This problem could be resolved, however, by integrating this equation along the characteristic 
u/p = const instead of the x-axis as we did to obtain eq.(16). The matter of bigger concern is 
that the increase of u x is obviously has to do with strong shock modification, so that P ~ 1. 
Clearly, under these circumstances the balance between the l.h.s. of eq.(14) and the pressure 
term on the r.h.s. leads to impossibly large I. Evidently, the second term on the r.h.s. must 
come into play before this has happened so that the steady state will be maintained by the 



14 


balance between this term and the pressure term while the l.h.s. will become sub-dominant. 
Thus, for / we have the following equation 

2 4i p - st ^ = ° < 28 > 

As it is often the case we may assume that the wave collision term St{I} oc I 2 and in 
the long wave limit k —» 0 it is also proportional to k 2 (that means to p~ 2 ). The pressure 
gradient may be estimated as P/L p where L p is the scale height of particles of momentum p 
which we assume for simplicity to be proportional to p as in the standard Bohrn case. Thus, 
for / we have 


1 2 ~ —pP (29) 

OL 14 

where a characterizes the strength of nonlinear wave interaction. Using the last estimate, 
instead of eq.(21) we have 


Tm dg_ 

x 

where L n \ = K 0 /uiy/aM A . Introducing the fast variable £ (22) with L n \ instead of L, and 
repeating the derivation in section 3.1, for g we obtain the following equation 



1 

3 


dg 

u x p—dx 

op 


(30) 


with the obvious solution 


w 3 + ^= t ;7 


ui dg 
y/gdZ 


S 


g — — tanh —-£ (31) 

w 2u\ 

This solution is valid for £ > > 0 (^o ~ G _1 ) whereas at — oo < £ < one should use the 

linear formula (16) for the wave spectral density and thus the solution (27) instead of eq.(31). 
The uniformly valid solution may be also obtained by using an interpolation between eq.(16) 
and eq.(29) for I. We will not need, however, the explicit form of the front transition in the 
particle distribution in the region £ ~ 1 which means x ~ Xf (p). We will merely exploit the 
fact that this transition is much narrower (its width is Ax ~ L n \/^G(xf,p)) than that of 
the main part G(x) in the interval Xf < x < oo. The spatial scale of the latter is at least 
~ kb/ui or even broader if the linear approximation (18) can be used, in which case the 
length scale is determined by the linear damping of Alfven waves (Drury et al. 1996). 
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The only characteristic of the above internal solution that is needed to calculate the 
external solution G(x,p) is the position of the front transition in g on the x,p- plane, i.e., 
we need to calculate the function x = Xf (p). To do this we return to eq.(20). We solved 
it by neglecting its l.h.s. and finally arrived at the result for g and thus for / in eq.(20) 
that contains the fast variable £. Generally, this produces large terms in the next order 
of approximation coming from the l.h.s. To avoid that we must choose the position of the 
transition front (£(x,p) = 0) in such a way that it coincides with one of the characteristics 
of the operator on the l.h.s. of eq.(20), i.e., 


d 


d 


u— + u x p— 
ox op 


z(x,p) = 0 


or 


Mp)-p^ = o 

The choice of the concrete characteristic is based on the existence of the absolute maximum 
momentum p max beyond which there are neither particles nor waves. That means 

P 

Uf(p ) = u [X{(p)] = Ui - 

Pmax 

Likewise, the function x = Xf (p) is defined as 


Xf (p) = u 



V 

Pmax 


3.3. External solution 

While having obtained the form and the position x = Xf (p) of the narrow front in the 
particle distribution g(x,p) we still need to calculate g to the right from the front where 
it decays with x. This would be the external solution G(x,p ) introduced in the previous 
subsections, ft is clear that 


max g(x,p) « G(xf,p) = G 0 (p) 

X 

so that from eq.(21) we have the following equation 


1 

3 



d G 



dx 


u,(p)G„(p) 


(32) 
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The most important information about G(x,p) is contained in Go(p) for which from the last 
equation we obtain 


^-v(p)G 0 (p) + 4— G 0 (p) = 0 

op pi 

where we have introduced v(p) by 


(33) 


1 f Ul 

y (p) = tttT / G{x,p)du(x) 
Go{P) J U f(p) 

Eq.(33) can be easily solved for G 0 


(34) 


Go ijp) 


C 

v(p) 


exp 


u 1 


Pn 



(35) 


where C is a normalization constant that should be determined from matching this solution 
with that in the region p < p*. However, the function v depends on the solution itself. 
Fortunately, this quantity can be calculated prior to determining Go and therefore, this 
solution may be written in a closed form. To illustrate this, let us consider a particularly 
simple case of p ~ p max , and we will turn to the general case afterwards. Clearly, p ~ Pmax 
means U{(p) — u\ . Evidently, we may replace G in eq.(34) by G 0 so that for z/(p) we have 
v(p) ~ Ui — «f(p) = Mi (1 — p/pmax)- Thus, from eq.(35) we obtain the following shape of 
the cut-off near p max 


Go — C (p max — p) 3 (36) 

In the rest of the x, p-domain where x > Xf (p) and p is not close to p max , we may assume that 
the CR diffusion coefficient is close to its Bohm value. Indeed, in contrast to the phase space 
region x « X{ (p) at each given x,p there are waves generated along the entire characteristic 
of eq.(14) passing through this point of the phase space and occupying an extended region 
of the CR precursor, Figure 2. We may use then the asymptotic high Mach number solution 
found in (Maikov 1997) 


g(x,p) = 9o(p) exp 


1 + f3 

n(p) 


udx 


(37) 


Here /3 is numerically small (typically ~ 1/6) and this solution without /3-term manifests the 
balance between the diffusion and convection terms on the l.h.s. of eq.(13) which is more 
accurate approximation far upstream where the flow modification (r.h.s.) is weak. The flow 
profile depends on the form of ft(p) and for k — Kb oc p in the internal part of the shock 
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transition u(x) behaves linearly with x. Adopting this solution to the region x > Xf, we may 
write 


so that for v we have 


G(x,p) = G 0 (p) exp 



v(p) 



du exp 


/ 1 + P f u u'du' \ 
V «B Jut u x{u')J 


In the Bohm case we can use the simplified linear approximation for u(x) Maikov (1997), 
u = Uq + U\x/L c where L c = 7tkb(p)/26 ) -u 1 , 6 ~ 1.09 and it is implied that the maximum 
contribution to the particle pressure comes from the momentum p = p (specified later). Now 
we can express v in the form of an error integral 


v (p) ~ / duex p 




(1 + (3) L c , 2 _ 2\ 
2 u lKB (p) 1 {) 


(38) 


The algebra further simplihes in two limiting cases (the second of which has been already 
mentioned) 


v(p) = U\ 


7TUift B /2 (1 + P) L c: 

(1 -p/Pmax) , 


P 'C Pmax 
P — Pmax 


This yields the following asymptotic behaviour of Go(p) 


G 0 (p) 



Vp ( 1 + p) /Opexp ( K ~^V( 1 + P)pp/ 6 ) r 

(Pmax p) 


P ^ Pmax 
P — Pmax 


(39) 


This result was obtained for particles with momenta p > p* = p max /R = p m ax u o/ u \, whereas 
for p < p* we may use the spectra tabulated in (Maikov & Drury 2001) for different p max ■> 
that should be associated now with p = p . The matching of these two spectra should give 
the normalization constant C in the solution (39). This will be the subject of the next 
section. 


3.4. Connection with the main part of the spectrum 

A typical solution of the nonlinear acceleration problem with a prescribed maximum 
momentum p max calculated using the method of integral equations developed in (Maikov 
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1997; Maikov & Drury 2001) is shown in Figure 4 with the dash-dotted line. Since the 
influence of shock modification on the injection rate is not known for high p max (see, however, 
Kang et al. 2001, where the values of p rnax ~ 10 have been reached) we have taken the 
injection rate v « 0.1, i.e., well inside the interval between the points A and C on Figure 1 
(see section 2.1). 

To calculate an integral spectrum containing both the part modified by the wave com¬ 
pression in the shock precursor as well as its lower energy downstream part (p < p*) we 
proceed as follows (see eq.(37) for the spatial structure of the spectrum). In the momentum 
range p < p* = p ma x/R we can obviously use the same method of integral equations. How¬ 
ever, the role of maximum momentum is now played not by p m ax but by p ( i.e., a dynamical 
cut-off where the maximum contribution to the CR pressure is coming from). Furthermore, 
given v and p, we calculate the self-consistent flow structure with the precompression R 
shown in Figure 1 as well as the particle spectrum. The latter is shown in Figure 4 with the 
dashed line. Note, that the spectrum matching point p* with the cut-off area p* < p < p max 
is now determined and we are ready to obtain the final spectrum by calculating its cut-off 
part from eqs.(35,39,38). The spectrum is drawn with the full line. The momentum p may 
be obtained now as a point of maximum of the function pG o (particle partial pressure per 
logarithmic interval) from eq.(39) (upper line) which yields 


P 



0.1 p 


max 


It should be clear from our treatment that this formula is valid only when the shock is 
strongly modified, namely whenp* = p max jR < 0.1 p max (the case we are interested in). Since 
R cannot exceed M 3 / 4 ; this means that the shock must be sufficiently strong, M 3 / 4 > 10. 


4. Discussion 

There are at least two reasons to believe that the standard acceleration theory may have 
estimated the maximum particle energy or the form of the spectrum below it incorrectly. 
The first reason is simply a possible conflict with the observations of TeV-emission from 
SNRs as we discussed in the Introduction section. The second reason is a theoretical one, 
that arises naturally from considering the nonlinear response of the shock structure to the 
acceleration which is exemplified in Figure 1. According to this picture, the response is so 
strong that it is unlikely that the acceleration can proceed at the same rate with no change in 
physics after such a dramatic shock restructuring (pre-compression R may rise by 1-2 orders 
of magnitudes depending on the Mach number). Time dependent numerical simulations 
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(e.g., Kang et al. 2002) show that the modifications occur very quickly, and compression is 
increased substantially even before p ma x ~ 1 (note that this would be consistent with the 
bifurcation diagram in Figure 1 for initial u ~ ( c/ui)ncii/ni > 0.1, where ncii/ni is the 
ratio of CR number density at the shock to that of the background plasma far upstream). 
The shock modification, in turn, must follow rather abruptly after the maximum momentum 
has passed through the critical value, ft was argued recently (Maikov et al. 2000) that 
this should drive crucial acceleration parameters such as the maximum momentum and 
injection rate back to their critical values which must limit shock modification and settle 
it at some marginal level, the so called self-organized critical (SOC) level (see also Jones 
2000; Maikov & Drury 2001; Maikov & Diamond 2001 for more discussions of the critical 
interrelation between the injection, maximum energy and shock structure). Mathematically, 
the SOC state is characterized by the requirement of merging of the two critical points on 
the bifurcation diagram in Figure 1 into one inflection point on the v(R) graph. Perhaps the 
most appealing aspect of this approach is its ability to predict the values of all three order 
parameters (injection rate, maximum momentum and compression ratio) given the only 
control parameter (the Mach number) just from our knowledge of the nonlinear response 
Riy^Pmax) shown in Figure 1, and no further calculations. 

However, the required backreaction mechanisms on the injection and maximum momen¬ 
tum need to be demonstrated to operate. We have already discussed at a qualitative level 
how the injection rate is reduced by shock modification. The subject of this paper has been 
the reduction of particle momenta related to the formation of a spectral break at p = p*, as a 
result of wave compression in a modified shock precursor. The position of the spectral break 
is universally related to the degree of system nonlinearity R, since p* = p m ax/R■ Hence, the 
problem seems to be converted to the study of nonlinear properties of the acceleration that 
are formally known from the analytic solution shown in Figure 1. However, the injection 
rate v that is now required for accurate determination of the spectral break p* through R , 
may currently be obtained only from the SOC ansatz. ft should be also mentioned that 
strong reduction of p* is obviously not to be expected in oblique shocks, where the resonance 
relation kp oc B is approximately preserved due to the compression of B simultaneously with 
k. 


An equally important problem is that strong losses of particles between p* and p max 
must slow down the growth of Pmaxif) due to the reduction of resonant waves. As we 
argued in section 2.1, this may result in an order of magnitude slower acceleration than one 
would expect from the standard Bohm diffusion paradigm. Consequently the dynamically 
and observationally significant spectral break p* may be at least two orders of magnitude 
below the maximum momentum p max (again, depending on M ) that could be reached in 
the unimpeded acceleration which is normally implied in estimates of maximum energy 
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achievable in SNRs over their active life time. 

In addition to the above mentioned uncertainty in p ma x(t)i its relation to the position 
of the spectral break p* = Pmax/R also needs further clarification. Indeed, since R depends 
on a dynamical cut-off p which, in general, is linked to p* and p m ax, the latter relation is still 
implicit. It can be easily resolved, however, in a supercritical regime (the saturated part of 
the R{v) dependence in Figure 1, see also Maikov & Drury 2001 for details), which requires 1 
vp/pinj 3> M 3 / 4 . One simply has then R « M 3 / 4 . As it was argued, however, the injection 
is unlikely to be high enough to reach this regime. An additional argument against it is that 
the spectral break becomes unrealistically small in the M —> oo limit, since p* = p max /M 3//4 . 
In the opposite case vp/p tnj <C M 3 / 4 , the compression rate saturates at R ~ vp/p in j. Note 
that the injection rate must be still above critical, otherwise R ~ 1. Now we need to specify 
p. The simple approximation used in the previous section yielded p ~ 0.1p max , so that 
p* ~ 10 Pinj/v (independent of p max ) which may be regarded as a lower bound on p*. Indeed, 
the above relation between p and p* may be applied only to the outermost part of the shock 
transition (see Secs. 3.3,3.4). Downstream, the spectrum cuts off very sharply immediately 
beyond p*, section 3.1. Therefore, the dynamical cut-off p ~ p* and we obtain the following 

upper bound on p*, p* « a J PinjPmaxit) / v . 

It should be clear that unless v is dramatically reduced as a result of shock modification, 
even this upper bound places p* way below p max ■ This may be the reason for non-detection 
of protons at TeV energies in SNRs. Finally, this does not contradict to the detection of 
10-100 TeV electrons in e.g., SNR 1006 since they may be accelerated by other mechanisms 
(e.g., Papadopoulos 1981; Galeev 1984; Bykov & Uvarov 1999; Laming 2001) or may have 
higher radiation efficiency. 


This study is supported through UCSD by the US Department of Energy, Grant No. 
FG03-88ER53275. At the LIniversity of Minnesota this work is supported by NASA through 
grant NAG5-8428, and by the LIniversity of Minnesota Supercomputing Institute. 


lr This is strictly valid for n(p) oc p. 
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Fig. 1.— Response of the shock structure (bifurcation diagram) to the injection of thermal 
particles at the rate u. The strength of the response is characterized by the precompression 
of the flow in the CR shock precursor R = U\/uq. The flow Mach number M = 150. Different 
curves correspond to different values of maximum momentum normalized to me. For each 
given v and p ma xi one (for p max < Per — 500) or three (for p m ax > Per) solutions exist. Note 
that solution multiplicity does not exist for shocks with M < M cr ~ 70 (Maikov et al. 2000; 
Maikov & Drury 2001). Given an initial injection v and compression R at point A (with 
R(A ) « 1), the injection and R at point C is calculated as v{C) = R(C)u(A) (see text for 
further explanations). 
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Fig. 3.— The phase plane of accelerated particles in the flow velocity-particle momentum 
coordinates. The particles are localized between the heavy line (above which there are 
no resonant waves to confine them) and the light line where the particle density decays 
exponentially towards higher u (see text). The relevant transport processes are indicated by 


arrows. 










Fig. 4.— Particle spectra at a strong shock obtained from analytic solution (Maikov 1997; 
Maikov & Drury 2001) for M = 150 (as in Figure 1) and the injection rate v « 0.1. The 
dash-dotted line shows the solution with the abrupt momentum cut-off at p — p max — 10 5 . 
The spectrum drawn with the full line demonstrates the effect of wave compression calculated 
using formulae (35) and (38). The spectrum that would be obtained using the same technique 
as for the dash-dotted case but with the maximum momentum at p = p (see text) is shown 
by the dashed line. 










